module projectors
  !
  ! Projection generator
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
#define green(i) (iand(mana,i)/=0)
  use suty
  use ground_control
  use PLANAR_PEM
  implicit none
  private
  public init_proj
  public init_subset
  public save_Cj
  public load_Cj
  public enable_subsets
  public disable_subsets
  public bproj_2d
  public fproj,bproj
  pointer(pfproj,fproj)
  pointer(pbproj,bproj)
  interface 
     subroutine fproj(a,b)
       real*4,intent(in)::a(:,:,:)
       real*4,intent(inout)::b(:,:,:,:)
     end subroutine fproj
  end interface
  interface 
     subroutine bproj(a,b)
       real*4,intent(in)::a(:,:,:,:)
       real*4,intent(inout)::b(:,:,:)
     end subroutine bproj
  end interface
  !
  integer,public,parameter::OPT_DOUBLE = int(z'0001000000000000',8)
  !
  integer opt
  integer nth 
  integer nbuf
  integer nsub
  integer*2,allocatable,target::jxs(:,:),jys(:,:),jzs(:,:)
  real*4,allocatable,target::ws(:,:)
  !
  integer ndj,ndjn
  real*4,allocatable,target::D_j_n(:,:,:,:,:)
  real*4,allocatable,target::D_j  (:,:,:,:)
  real*4,pointer::D_j_0(:,:,:,:)
  real*4,pointer::D_j_1(:,:,:,:)
  real*4,pointer::D_j_2(:,:,:,:)
  real*4,pointer::D_j_3(:,:,:,:)
  !
  integer npin
  real*4,allocatable,target::p_i_n(:,:,:,:,:)
  real*4,pointer::p_i(:,:,:,:) ! <<<
  real*4,pointer::p_i_0(:,:,:,:)
  real*4,pointer::p_i_1(:,:,:,:)
  real*4,pointer::p_i_2(:,:,:,:)
  real*4,pointer::p_i_3(:,:,:,:)
  !
  real*4,allocatable,target::q_i_n(:,:,:,:,:)
  real*4,pointer::q_i(:,:,:,:)
  real*4,pointer::q_i_0(:,:,:,:)
  real*4,pointer::q_i_1(:,:,:,:)
  real*4,pointer::q_i_2(:,:,:,:)
  real*4,pointer::q_i_3(:,:,:,:)
  !
  real*4,allocatable,target::x_j_n(:,:,:,:)
  real*4,pointer::x_j(:,:,:) ! <<<
  real*4,pointer::x_j_0(:,:,:)
  real*4,pointer::x_j_1(:,:,:)
  real*4,pointer::x_j_2(:,:,:)
  real*4,pointer::x_j_3(:,:,:)
  !
  integer,allocatable,target::offsets(:,:)
  integer,pointer::xoff(:),yoff(:)
  integer ox,oy,ux,uy,isub
  !
  type proj_arg
     integer n,no,nsx,nsy
     integer ix1,iy1,dx,dy
  end type proj_arg
  type(proj_arg),allocatable,target::pas(:)
  !
  external setma_d,setma_dr
  pointer(psetma_d,setma_d)
  pointer(psetma_dr,setma_dr)
  !
  external sumc,sumcr,sumx,sumxr,sumy,sumyr,sumxy,sumxyr
  pointer(psumc,sumc)
  pointer(psumcr,sumcr)
  pointer(psumx,sumx)
  pointer(psumxr,sumxr)
  pointer(psumy,sumy)
  pointer(psumyr,sumyr)
  pointer(psumxy,sumxy)
  pointer(psumxyr,sumxyr)
  !
contains
  !
  subroutine echo(m)
    character*(*),intent(in)::m
    if(.not.green(OPT_VERBOSE)) return
    write(*,10) trim(m)
    call flush(6)
10  format(x,"(pj) ",a)
  end subroutine echo
  !
  subroutine echo_red(mess,code)
    character*(*),intent(in)::mess
    integer,intent(in),optional::code
    integer ec
    character*(*),parameter::swan_song="*** emergency escape from projectors."
    write(0,*) trim(adjustl(mess))
    call set_rstat(RSTAT_ERROR)
    ec=-1
    if(present(code)) ec=code
    call exit(ec)
  end subroutine echo_red
  !
  subroutine echo_si(s,i)
    character*(*),intent(in)::s
    integer,intent(in)::i
    if(green(OPT_DEBUG)) &
         call echo(trim(s)//repeat(" ",8-len_trim(s))//"launched on thread "//trim(i2a(i)))
  end subroutine echo_si
  !
  subroutine echo_sx(t0)
    use omp_lib
    real,intent(in)::t0
    real t1
    if(green(OPT_DEBUG)) then
       call cpu_time(t1)
       call echo("Thread "//trim(i2a(omp_get_thread_num()+1))//" CPU time: "//trim(f2a(t1-t0,3)))
    end if
  end subroutine echo_sx
  !
  subroutine enable_subsets()
    if(nsub==1) return
    psumc   = loc(sumcs)
    psumcr  = loc(sumcrs)
    psumx   = loc(sumxs)
    psumxr  = loc(sumxrs)
    psumy   = loc(sumys)
    psumyr  = loc(sumyrs)
    psumxy  = loc(sumxys)
    psumxyr = loc(sumxyrs)
  end subroutine enable_subsets
  !
  subroutine disable_subsets()
    psumc   = loc(sumca)
    psumcr  = loc(sumcra)
    psumx   = loc(sumxa)
    psumxr  = loc(sumxra)
    psumy   = loc(sumya)
    psumyr  = loc(sumyra)
    psumxy  = loc(sumxya)
    psumxyr = loc(sumxyra)
  end subroutine disable_subsets
  !
  subroutine init_proj(ires,&
       num_buffer,num_threads,num_subsets,option)
    use cuty
    integer,intent(in)::ires
    integer,intent(in)::num_buffer
    integer,intent(in)::num_threads
    integer,intent(in)::num_subsets
    integer,intent(in),optional::option
    integer nnz
    opt=0
    !
    if(present(option)) opt=option
    nsub=num_subsets
    nbuf=num_buffer
    nth=num_threads
    !
    call disable_subsets
    !
    select case(ires)
    case(2)
       pfproj=loc(proj_sysmat_j)
       pbproj=loc(proj_sysmat_i)
    case(1)
       pfproj=loc(proj_sysmat_j2)
       pbproj=loc(proj_sysmat_i2)
    end select
    !
    call initenvc(nx,ny_true)
    call alloc_js_aligned(nbuf,nth)
    !
    if(allocated(D_j_n)) deallocate(D_j_n)
    nnz=nz+(nz-8*(nz/8)) ! <<<<<<<<<<<<<<<<<<
    allocate(D_j_n(nh,ny,nnz,nth,4)) ! <<<
    ndjn=nh*ny*(nnz)*nth*4 ! <<<
    D_j_0=>D_j_n(:,:,:,:,1)
    D_j_1=>D_j_n(:,:,:,:,2)
    D_j_2=>D_j_n(:,:,:,:,3)
    D_j_3=>D_j_n(:,:,:,:,4)
    if(allocated(D_j)) deallocate(D_j)
    allocate(D_j(nx,ny,nnz,nth)) ! <<<
    ndj=nx*ny*(nnz)*nth ! <<<
    !
    if(allocated(x_j_n)) deallocate(x_j_n)
    allocate(x_j_n(nh,ny,nnz,4)) ! <<<
    x_j_0=>x_j_n(:,:,:,1)
    x_j_1=>x_j_n(:,:,:,2)
    x_j_2=>x_j_n(:,:,:,3)
    x_j_3=>x_j_n(:,:,:,4)
    !
    if(allocated(p_i_n)) deallocate(p_i_n)
    allocate(p_i_n(mx,my,mx,my,2))
    p_i_0=>p_i_n(:,:,:,:,1)
    p_i_1=>p_i_n(:,:,:,:,1)
    p_i_2=>p_i_n(:,:,:,:,2)
    p_i_3=>p_i_n(:,:,:,:,2)
    npin=mxy2*2
    !
    if(allocated(q_i_n)) deallocate(q_i_n)
    allocate(q_i_n(mx,my,mx,my,2))
    q_i_0=>q_i_n(:,:,:,:,1)
    q_i_1=>q_i_n(:,:,:,:,1)
    q_i_2=>q_i_n(:,:,:,:,2)
    q_i_3=>q_i_n(:,:,:,:,2)
    !
    if(allocated(pas)) deallocate(pas)
    allocate(pas(nth))
    !
    select case(nsub)
    case(4)
       if(allocated(offsets)) deallocate(offsets)
       allocate(offsets(nsub,2))
       offsets=0
       xoff=>offsets(:,1)
       yoff=>offsets(:,2)
       xoff=[0,1,0,1]
       yoff=[0,0,1,1]
    case(1)       
    case default
       call echo_red("*** Error: unsupported number of subsets.")
    end select
    !
  end subroutine init_proj
  !
  subroutine alloc_js_aligned(n,m)
    use cuty
    integer,intent(in)::n,m
    integer n8,n4
    if(allocated(jxs)) deallocate(jxs)
    if(allocated(jys)) deallocate(jys)
    if(allocated(jzs)) deallocate(jzs)
    if(allocated(ws)) deallocate(ws)   
    if(m==1) then
       n8=n
       n4=n
       allocate(jxs(n,m),jys(n,m),jzs(n,m),ws(n,m))
    else
       n8=8*(n/8+1)
       n4=4*(n/4+1)
       allocate(jxs(n8,m),jys(n8,m),jzs(n8,m),ws(n4,m))
    end if
    call initenvd(n8,&
         loc(jxs),loc(jys),loc(jzs),loc(ws),n4)
  end subroutine alloc_js_aligned
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine bproj_2d(pi,xj,d)
    real*4,intent(in)::pi(:,:,:,:)
    real*4,intent(out),target::xj(:,:,:)
    real,intent(in)::d
    real jxy0(2),z1
    integer iy2_
    xj=0.0_4
    jxy0(1)=jx0
    jxy0(2)=jy0
    z1=-d/2.0_dp
    !$omp parallel do shared(xj) schedule(dynamic)
    do iy2_=2,my-1
       call do_proj(iy2_)
    end do
    !$omp end parallel do
  contains
    subroutine do_proj(iy2)
      integer,intent(in)::iy2
      real xy1(2),xy2(2),xy0(2),p,dxy(2)
      integer ix1,iy1,ix2,jxy(2),jz
      real f(2)
      xy2(2)=(2*iy2-my)-1
      do ix2=2,mx-1
         xy2(1)=(2*ix2-mx)-1
         do iy1=2,my-1
            xy1(2)=(2*iy1-my)-1
            do ix1=2,mx-1
               p=pi(ix1,iy1,ix2,iy2)
               if(p==0) cycle
               xy1(1)=(2*ix1-mx)-1
               dxy=xy2-xy1     
               do jz=1,nz
                  xy0=xy1+((jz-jz0)*dz-z1)/d*dxy
                  f=xy0+jxy0-int(xy0+jxy0)
                  jxy=int(xy0+jxy0)
                  if(jxy(1)<=0.or.jxy(2)<=0&
                       .or.jxy(1)>nx.or.jxy(2)>ny) cycle
                  xj(jxy(1),   jxy(2),jz)   =  xj(jxy(1),   jxy(2),jz)  +real(p*(1.0_dp-f(1))*(1.0_dp-f(2)),sp)
                  xj(jxy(1)+1, jxy(2),jz)   =  xj(jxy(1)+1, jxy(2),jz)  +real(p*f(1)*(1.0_dp-f(2)),sp)
                  xj(jxy(1),   jxy(2)+1,jz) =  xj(jxy(1),   jxy(2)+1,jz)+real(p*(1.0_dp-f(1))*f(2),sp)
                  xj(jxy(1)+1, jxy(2)+1,jz) =  xj(jxy(1)+1, jxy(2)+1,jz)+real(p*f(1)*f(2),sp)
               end do
            end do
         end do
      end do
    end subroutine do_proj
  end subroutine bproj_2d
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine load_Cj(f,cj)
    character*(*),intent(in)::f
    real*4,intent(out)::Cj(:,:,:)
    integer istat
    open(unit=1,file=f,iostat=istat,status="old",access="stream")
    if(istat/=0) call echo_red(mess_ioerr(EID_OPEN,f,istat))
    read(1,iostat=istat) Cj
    close(1)
    if(istat/=0) call echo_red(mess_ioerr(EID_READ,f,istat))    
    call echo("Loaded Cj: "//trim(f))
  end subroutine load_Cj
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine save_Cj(f,v)
    character*(*),intent(in)::f
    real*4,intent(in)::v(:,:,:)
    character(MAXPATHLEN) fout
    integer istat
    fout=trim(f)//"_"//trim(i2a(size(v,1)))//"x"&
         //trim(i2a(size(v,2)))//"x"&
         //trim(i2a(size(v,3)))//".cj"
    call write_v(fout,-1,-1,loc(v),size(v),istat,OFMT_RAW)
    if(istat/=0) call echo_red(mess_ioerr(EID_WRITE,fout,istat))
  end subroutine save_Cj
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine sumca
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMC",i)
    pa=>pas(i)
    pa%ix1=1; pa%iy1=1
    do dy2=0,mdxy_max-1
       iy2=1+dy2
       pa%nsy=my_true-dy2
       pa%dy=iy2
       do dx2=dy2,mdxy_max-1
          ix2=1+dx2
          sm=>sysmat(ix2,iy2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=ix2
          pa%no=sm%noy
          pa%nsx=mx-dx2
          call cpptra(pa%n,sm%q,i)
          call setma_dr(i)
       end do
    end do
  end subroutine sumca
  !
  subroutine sumcra
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMCr",i)
    pa=>pas(i)
    pa%ix1=1; pa%iy1=1
    do dy2=1,mdxy_max-1
       iy2=1+dy2
       pa%dy=iy2
       pa%nsy=my_true-dy2
       do dx2=0,dy2-1
          ix2=1+dx2
          sm=>sysmat(iy2,ix2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=ix2
          pa%no=sm%noy
          pa%nsx=mx-dx2
          call cpptrb(pa%n,sm%q,i)
          call setma_d(i)
       end do
    end do
  end subroutine sumcra
  !
  subroutine sumxa
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMx",i)
    pa=>pas(i)
    pa%ix1=mx; pa%iy1=1
    do dy2=0,mdxy_max-1
       iy2=1+dy2
       pa%dy=iy2
       pa%nsy=my_true-dy2
       do dx2=max(1,dy2),mdxy_max-1
          ix2=1+dx2
          sm=>sysmat(ix2,iy2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=mx-dx2+1
          pa%no=sm%noy
          pa%nsx=-mx+dx2
          call cpptrax(pa%n,sm%q,i)
          call setma_dr(i)
       end do
    end do
  end subroutine sumxa
  !
  subroutine sumxra
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMxr",i)
    pa=>pas(i)
    pa%ix1=mx; pa%iy1=1
    do dy2=2,mdxy_max-1
       iy2=1+dy2
       pa%dy=iy2
       pa%nsy=my_true-dy2
       do dx2=1,dy2-1
          ix2=1+dx2
          sm=>sysmat(iy2,ix2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=mx-dx2+1
          pa%no=sm%noy
          pa%nsx=-mx+dx2
          call cpptrbx(pa%n,sm%q,i)
          call setma_d(i)
       end do
    end do
  end subroutine sumxra
  !
  subroutine sumya
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMy",i)
    pa=>pas(i)
    pa%ix1=1; pa%iy1=my_true
    do dy2=1,mdxy_max-1
       iy2=1+dy2
       pa%dy=my-dy2+1
       pa%nsy=-my_true+dy2
       do dx2=dy2,mdxy_max-1
          ix2=1+dx2
          sm=>sysmat(ix2,iy2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=ix2
          pa%no=sm%noy
          pa%nsx=mx-dx2
          call cpptray(pa%n,sm%q,i)
          call setma_dr(i)
       end do
    end do
  end subroutine sumya
  !
  subroutine sumyra
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMyr",i)
    pa=>pas(i)
    pa%ix1=1; pa%iy1=my_true
    do dy2=1,mdxy_max-1
       iy2=1+dy2
       pa%dy=my-dy2+1
       pa%nsy=-my_true+dy2
       do dx2=0,dy2-1
          ix2=1+dx2
          sm=>sysmat(iy2,ix2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=ix2
          pa%no=sm%noy
          pa%nsx=mx-dx2
          call cpptrby(pa%n,sm%q,i)
          call setma_d(i)
       end do
    end do
  end subroutine sumyra
  !
  subroutine sumxya
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMxy",i)
    pa=>pas(i)
    pa%ix1=mx; pa%iy1=my_true
    do dy2=1,mdxy_max-1
       iy2=1+dy2
       pa%dy=my-dy2+1
       pa%nsy=-my_true+dy2
       do dx2=dy2,mdxy_max-1
          ix2=1+dx2
          sm=>sysmat(ix2,iy2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=mx-dx2+1
          pa%no=sm%noy
          pa%nsx=-mx+dx2
          call cpptraxy(pa%n,sm%q,i)
          call setma_dr(i)
       end do
    end do
  end subroutine sumxya
  !
  subroutine sumxyra
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMxyr",i)
    pa=>pas(i)
    pa%ix1=mx; pa%iy1=my_true
    do dy2=2,mdxy_max-1
       iy2=1+dy2
       pa%dy=my-dy2+1
       pa%nsy=-my_true+dy2
       do dx2=1,dy2-1
          ix2=1+dx2
          sm=>sysmat(iy2,ix2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=mx-dx2+1
          pa%no=sm%noy
          pa%nsx=-mx+dx2
          call cpptrbxy(pa%n,sm%q,i)
          call setma_d(i)
       end do
    end do
  end subroutine sumxyra
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine run_sum
    real t0,t1 
    call cpu_time(t0)
    !$omp parallel
    !$omp sections 
    call sumc
    call echo_sx(t0)
    !
    !$omp section
    call sumcr
    call echo_sx(t0)
    !
    !$omp section
    call sumx
    call echo_sx(t0)
    !
    !$omp section
    call sumxr
    call echo_sx(t0)
    !
    !$omp section
    call sumy
    call echo_sx(t1-t0)
    !
    !$omp section
    call sumyr
    call echo_sx(t0)
    !
    !$omp section
    call sumxy
    call echo_sx(t0)
    !
    !$omp section
    call sumxyr
    call echo_sx(t0)
    !$omp end sections
    !$omp end parallel
  end subroutine run_sum
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine setpi_d(i)
    use vuty
    integer,intent(in)::i
    integer k
    integer po,pe
    real s
    type(proj_arg),pointer::pa
    pa=>pas(i)
    po=loc(p_i_1(pa%ix1,pa%iy1,pa%dx,pa%dy))
    pe=loc(p_i_0(pa%ix1,pa%iy1,pa%dx,pa%dy))
    !$omp parallel private(k)
    !$omp sections
    do k=1,pa%no
       call muladdpsl(pa%nsx,loc(x_j_1((jxs(k,i)+1)/2,jys(k,i),jzs(k,i))),&
            loc(ws(k,i)),pa%nsy,po)
    end do
    !$omp section
    do k=pa%no+1,pa%n
       call muladdpsl(pa%nsx,loc(x_j_0(jxs(k,i)/2,jys(k,i),jzs(k,i))),&
            loc(ws(k,i)),pa%nsy,pe)
    end do
    !$omp section
    s=0
    do k=1,pa%n
       s=s+ws(k,i)*x_j(jxs(k,i),jys(k,i),jzs(k,i))
    end do
    p_i_0(pa%ix1,pa%iy1,pa%dx,pa%dy)=real(s,4)
    !$omp end sections
    !$omp end parallel
  end subroutine setpi_d
  !
  subroutine setpi_dr(i)
    use vuty
    integer,intent(in)::i
    integer k
    integer po,pe
    real s
    type(proj_arg),pointer::pa
    pa=>pas(i)
    po=loc(p_i_3(pa%iy1,pa%ix1,pa%dy,pa%dx))
    pe=loc(p_i_2(pa%iy1,pa%ix1,pa%dy,pa%dx))
    !$omp parallel private(k)
    !$omp sections
    do k=1,pa%no
       call muladdpsl(pa%nsy,loc(x_j_3((jys(k,i)+1)/2,jxs(k,i),jzs(k,i))),&
            loc(ws(k,i)),pa%nsx,po)
    end do
    !$omp section
    do k=pa%no+1,pa%n
       call muladdpsl(pa%nsy,loc(x_j_2(jys(k,i)/2,jxs(k,i),jzs(k,i))),&
            loc(ws(k,i)),pa%nsx,pe)
    end do
    !$omp section
    s=0
    do k=1,pa%n
       s=s+ws(k,i)*x_j(jxs(k,i),jys(k,i),jzs(k,i))
    end do
    p_i_0(pa%ix1,pa%iy1,pa%dx,pa%dy)=real(s,4)
    !$omp end sections
    !$omp end parallel
  end subroutine setpi_dr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine proj_sysmat_j(xj_in,pi)
    use omp_lib
    real*4,intent(out),target::pi(:,:,:,:)
    real*4,intent(inout),target::xj_in(:,:,:)
    real*4,pointer::xj(:,:,:)
    real*4,target::tmpimg(nx,ny,nz)
    real t0,t1
    integer rt(8,2)
    !
    call echo("proj_sysmat_j: begin")
    call date_and_time(values=rt(:,1))
    !
    call cpu_time(t0)
    !
    if(green(OPT_USE_FILTER)) then
       tmpimg=xj_in
       call conv_img(tmpimg,.FALSE.)
       xj=>tmpimg
    else
       xj=>xj_in
    end if
    !
    call cle2d(loc(p_i_n),npin)
    call sepxj
    x_j=>xj
    psetma_d =loc(setpi_d)
    psetma_dr=loc(setpi_dr)
    call cpu_time(t1)
    call echo_sx(t1-t0)
    !
    call run_sum
    !
    call cpu_time(t0)
    call mergepi
    call cpu_time(t1)
    call echo_sx(t1-t0)
    !
    call date_and_time(values=rt(:,2))
    call echo("proj_sysmat_j: end: "//trim(f2a(etime(rt(:,1),rt(:,2)),3)))
    !
  contains
    !
    subroutine sepxj
      integer x,y,z
      integer n,xodd,yodd
      xodd=iand(nx,1)
      yodd=iand(ny,1)
      n=nh
      if(xodd==1) n=nh-1
      do z=1,nz
         do y=1,ny
            do x=1,n
               x_j_1(x,y,z)=xj(2*x-1,y,z)
               x_j_0(x,y,z)=xj(2*x,  y,z)
            end do
            if(xodd==1) x_j_1(nh,y,z)=xj(nx,y,z)
         end do
      end do
      do z=1,nz
         do x=1,nx
            do y=1,n
               x_j_3(y,x,z)=xj(x,2*y-1,z)
               x_j_2(y,x,z)=xj(x,  2*y,z)
            end do
            if(yodd==1) x_j_3(nh,x,z)=xj(x,ny,z)
         end do
      end do
    end subroutine sepxj
    !
    subroutine mergepi
      integer x,y,ddx,ddy
      !$omp parallel private(x,y,ddx,ddy)
      !$omp sections
      do ddy=0,mdxy_max-1
         do ddx=0,mdxy_max-1 ! positive dx
            do y=1,my-ddy
               do x=1,mx-ddx
                  pi(x,y,ddx+1,ddy+1)=&
                       p_i_0(x,y,ddx+1,ddy+1)+p_i_2(y,x,ddy+1,ddx+1)
               end do
            end do
         end do
         do ddx=1,mdxy_max ! negative dx
            do y=1,my-ddy
               do x=ddx+1,mx
                  pi(x,y,mx-ddx+1,ddy+1)=&
                       p_i_0(x,y,mx-ddx+1,ddy+1)&
                       +p_i_2(y,x,ddy+1,mx-ddx+1)
               end do
            end do
         end do
      end do
      !$omp section
      do ddy=1,mdxy_max-1 ! negative dy
         do ddx=0,mdxy_max-1 ! positive dx
            do y=ddy+1,my
               do x=1,mx-ddx
                  pi(x,y,ddx+1,my-ddy+1)=&
                       p_i_0(x,y,ddx+1,my-ddy+1)&
                       +p_i_2(y,x,my-ddy+1,ddx+1)
               end do
            end do
         end do
         do ddx=1,mdxy_max-1 ! negative dx
            do y=ddy+1,my
               do x=ddx+1,mx
                  pi(x,y,mx-ddx+1,my-ddy+1)=&
                       p_i_0(x,y,my-ddx+1,my-ddy+1)&
                       +p_i_2(y,x,my-ddy+1,my-ddx+1)
               end do
            end do
         end do
      end do
      !$omp end sections
      !$omp end parallel
    end subroutine mergepi
    !
  end subroutine proj_sysmat_j
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine setpi_d2(i)
    use vuty
    integer,intent(in)::i
    integer k
    integer p
    real s
    type(proj_arg),pointer::pa
    pa=>pas(i)
    p=loc(p_i(pa%ix1,pa%iy1,pa%dx,pa%dy))
    do k=1,pa%n
       call muladdpsl2(pa%nsx,loc(x_j(jxs(k,i),jys(k,i),jzs(k,i))),&
            loc(ws(k,i)),pa%nsy,p)
    end do
    s=0
    do k=1,pa%n
       s=s+ws(k,i)*x_j(jxs(k,i),jys(k,i),jzs(k,i))
    end do
    p_i(pa%ix1,pa%iy1,pa%dx,pa%dy)=real(s,4)
  end subroutine setpi_d2
  !
  subroutine setpi_dr2(i)
    use vuty
    integer,intent(in)::i
    integer k
    integer p
    real s
    type(proj_arg),pointer::pa
    pa=>pas(i)
    p=loc(p_i(pa%iy1,pa%ix1,pa%dy,pa%dx))
    do k=1,pa%n
       call muladdpsl2(pa%nsy,loc(x_j(jys(k,i),jxs(k,i),jzs(k,i))),&
            loc(ws(k,i)),pa%nsx,p)
    end do
    s=0
    do k=1,pa%n
       s=s+ws(k,i)*x_j(jxs(k,i),jys(k,i),jzs(k,i))
    end do
    p_i(pa%ix1,pa%iy1,pa%dx,pa%dy)=real(s,4)
  end subroutine setpi_dr2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine proj_sysmat_j2(xj,pi)
    real*4,intent(out),target::pi(:,:,:,:)
    real*4,intent(inout),target::xj(:,:,:)
    real t0,t1
    integer rt(8,2)
    !
    call echo("proj_sysmat_j2: begin")
    call date_and_time(values=rt(:,1))
    !
    call cpu_time(t0)
    call cle2d(loc(pi),size(pi))
    x_j=>xj
    p_i=>pi
    psetma_d =loc(setpi_d2)
    psetma_dr=loc(setpi_dr2)
    call cpu_time(t1)
    call echo_sx(t1-t0)
    !
    call run_sum
    !
    call date_and_time(values=rt(:,2))
    call echo("proj_sysmat_j2: end: "//trim(f2a(etime(rt(:,1),rt(:,2)),3)))
  end subroutine proj_sysmat_j2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine setdj_d(i)
    use vuty
    integer,intent(in)::i
    integer k
    integer po,pe
    real*4,pointer::D_j_o(:,:,:),D_j_e(:,:,:),D_j_i(:,:,:)
    real q
    type(proj_arg),pointer::pa
    pa=>pas(i)
    po=loc(q_i_1(pa%ix1,pa%iy1,pa%dx,pa%dy))
    pe=loc(q_i_0(pa%ix1,pa%iy1,pa%dx,pa%dy))
    D_j_o=>D_j_1(:,:,:,i)
    D_j_e=>D_j_0(:,:,:,i)
    D_j_i=>D_j(:,:,:,i)
    !$omp parallel private(k)
    !$omp sections
    do k=1,pa%no
       call muladdps(pa%nsx,po,loc(ws(k,i)),pa%nsy,&
            loc(D_j_o((jxs(k,i)+1)/2,jys(k,i),jzs(k,i))))
    end do
    !$omp section
    do k=pa%no+1,pa%n
       call muladdps(pa%nsx,pe,loc(ws(k,i)),pa%nsy,&
            loc(D_j_e(jxs(k,i)/2,jys(k,i),jzs(k,i))))
    end do
    !$omp section
    q=q_i_0(pa%ix1,pa%iy1,pa%dx,pa%dy)
    do k=1,pa%n
       D_j_i(jxs(k,i),jys(k,i),jzs(k,i))=&
            D_j_i(jxs(k,i),jys(k,i),jzs(k,i))+real(ws(k,i)*q,4)
    end do
    !$omp end sections
    !$omp end parallel
  end subroutine setdj_d
  !
  subroutine setdj_dr(i)
    use vuty
    integer,intent(in)::i
    integer k
    integer po,pe
    real*4,pointer::D_j_o(:,:,:),D_j_e(:,:,:),D_j_i(:,:,:)
    real q
    type(proj_arg),pointer::pa
    pa=>pas(i)
    po=loc(q_i_3(pa%iy1,pa%ix1,pa%dy,pa%dx))
    pe=loc(q_i_2(pa%iy1,pa%ix1,pa%dy,pa%dx))
    D_j_o=>D_j_3(:,:,:,i)
    D_j_e=>D_j_2(:,:,:,i)
    D_j_i=>D_j(:,:,:,i)
    !$omp parallel private(k)
    !$omp sections
    do k=1,pa%no
       call muladdps(pa%nsy,po,loc(ws(k,i)),pa%nsx,&
            loc(D_j_o((jys(k,i)+1)/2,jxs(k,i),jzs(k,i))))
    end do
    !$omp section
    do k=pa%no+1,pa%n
       call muladdps(pa%nsy,pe,loc(ws(k,i)),pa%nsx,&
            loc(D_j_e(jys(k,i)/2,jxs(k,i),jzs(k,i))))
    end do
    !$omp section
    q=q_i_0(pa%ix1,pa%iy1,pa%dx,pa%dy)
    do k=1,pa%n
       D_j_i(jxs(k,i),jys(k,i),jzs(k,i))=&
            D_j_i(jxs(k,i),jys(k,i),jzs(k,i))+real(ws(k,i)*q,4)
    end do
    !$omp end sections
    !$omp end parallel
  end subroutine setdj_dr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine proj_sysmat_i(qi,Dj)
    real*4,intent(inout)::qi(:,:,:,:)
    real*4,intent(out)::Dj(:,:,:)
    real t0,t1
    integer rt(8,2)
    !
    call echo("proj_sysmat_i: begin")
    call date_and_time(values=rt(:,1))
    !
    call cpu_time(t0)
    call cle2d(loc(D_j_n),ndjn)
    call cle2d(loc(D_j),ndj)
    call sepqi
    psetma_d =loc(setdj_d)
    psetma_dr=loc(setdj_dr)
    call cpu_time(t1)
    call echo_sx(t1-t0)
    !
    call run_sum

    !
    call cpu_time(t0)
    call mergedj

    if(green(OPT_USE_FILTER)) call conv_img(Dj,.false.)

    call cpu_time(t1)
    call echo_sx(t1-t0)
    !
    call date_and_time(values=rt(:,2))
    call echo("proj_sysmat_i: end: "//trim(f2a(etime(rt(:,1),rt(:,2)),3)))
    !
  contains
    !
    subroutine mergedj
      integer jjx,jjy,jjz
      integer xodd,yodd
      integer n
!      forall(jjx=1:nh,jjy=1:ny,jjz=1:nz)  ! FIX20090324
      forall(jjx=1:nh,jjy=1:ny,jjz=1:nz) ! BUG20090324
         D_j_0(jjx,jjy,jjz,1)=sum(D_j_0(jjx,jjy,jjz,:))
         D_j_1(jjx,jjy,jjz,1)=sum(D_j_1(jjx,jjy,jjz,:))
         D_j_2(jjx,jjy,jjz,1)=sum(D_j_2(jjx,jjy,jjz,:))
         D_j_3(jjx,jjy,jjz,1)=sum(D_j_3(jjx,jjy,jjz,:))
      end forall
      forall(jjx=1:nx,jjy=1:ny,jjz=1:nz) ! BUG20090324
         D_j(jjx,jjy,jjz,1)  =sum(D_j(jjx,jjy,jjz,:))
      end forall
      call cle2d(loc(Dj),nv)
      n=nh
      xodd=iand(nx,1)
      yodd=iand(ny,1)
      if(xodd==1) n=n-1
      do jjz=1,nz
         do jjy=1,ny
            do jjx=1,n
               Dj(2*jjx-1,jjy,jjz) = D_j_1(jjx,jjy,jjz,1)
               Dj(2*jjx,  jjy,jjz) = D_j_0(jjx,jjy,jjz,1)
            end do
            if(xodd==1) Dj(nx,jjy,jjz) = D_j_1(nh,jjy,jjz,1)
         end do
      end do
      do jjz=1,nz
         do jjx=1,nx
            do jjy=1,n
               Dj(jjx,2*jjy-1,jjz) = Dj(jjx,2*jjy-1,jjz)+D_j_3(jjy,jjx,jjz,1)
               Dj(jjx,  2*jjy,jjz) = Dj(jjx,  2*jjy,jjz)+D_j_2(jjy,jjx,jjz,1)
            end do
            if(yodd==1) Dj(jjx,ny,jjz) = Dj(jjx,ny,jjz)+D_j_3(nh,jjx,jjz,1)
         end do
      end do
      forall(jjx=1:nx,jjy=1:ny,jjz=1:nz)
         Dj(jjx,jjy,jjz)=Dj(jjx,jjy,jjz)+D_j(jjx,jjy,jjz,1)
      end forall
    end subroutine mergedj
    !
    subroutine sepqi
      integer x,y,ddx,ddy
      call cp2d(loc(q_i_0),loc(qi),mxy2)
      do ddy=1,my
         do ddx=1,mx
            do y=1,my
               do x=1,mx
                  q_i_2(x,y,ddx,ddy)=qi(y,x,ddy,ddx)
               end do
            end do
         end do
      end do
    end subroutine sepqi
    !
  end subroutine proj_sysmat_i
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine setdj_d2(i)
    use vuty
    integer,intent(in)::i
    integer k
    integer p
    real*4,pointer::D_j_i(:,:,:)
    real q
    type(proj_arg),pointer::pa
    pa=>pas(i)
    p=loc(q_i(pa%ix1,pa%iy1,pa%dx,pa%dy))
    D_j_i=>D_j(:,:,:,i)
    do k=1,pa%n
       call muladdps2(pa%nsx,p,loc(ws(k,i)),pa%nsy,&
            loc(D_j_i(jxs(k,i),jys(k,i),jzs(k,i))))
    end do
    q=q_i(pa%ix1,pa%iy1,pa%dx,pa%dy)
    do k=1,pa%n
       D_j_i(jxs(k,i),jys(k,i),jzs(k,i))=&
            D_j_i(jxs(k,i),jys(k,i),jzs(k,i))+real(ws(k,i)*q,4)
    end do
  end subroutine setdj_d2
  !
  subroutine setdj_dr2(i)
    use vuty
    integer,intent(in)::i
    integer k
    integer p
    real*4,pointer::D_j_i(:,:,:)
    real q
    type(proj_arg),pointer::pa
    pa=>pas(i)
    p=loc(q_i(pa%iy1,pa%ix1,pa%dy,pa%dx))
    D_j_i=>D_j(:,:,:,i)
    do k=1,pa%n
       call muladdps2(pa%nsy,p,loc(ws(k,i)),pa%nsx,&
            loc(D_j_i(jys(k,i),jxs(k,i),jzs(k,i))))
    end do
    q=q_i(pa%ix1,pa%iy1,pa%dx,pa%dy)
    do k=1,pa%n
       D_j_i(jxs(k,i),jys(k,i),jzs(k,i))=&
            D_j_i(jxs(k,i),jys(k,i),jzs(k,i))+real(ws(k,i)*q,4)
    end do
  end subroutine setdj_dr2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine proj_sysmat_i2(qi,Dj)
    real*4,intent(inout),target::qi(:,:,:,:)
    real*4,intent(out)::Dj(:,:,:)
    real t0,t1
    integer rt(8,2)
    !
    call echo("proj_sysmat_i2: begin")
    call date_and_time(values=rt(:,1))
    !
    call cpu_time(t0)
    call cle2d(loc(D_j),ndj)
    q_i=>qi
    psetma_d =loc(setdj_d2)
    psetma_dr=loc(setdj_dr2)
    call cpu_time(t1)
    call echo_sx(t1-t0)
    !
    call run_sum
    call mergedj
    !    
    call date_and_time(values=rt(:,2))
    call echo("proj_sysmat_i2: end: "//trim(f2a(etime(rt(:,1),rt(:,2)),3)))
    !
  contains
    subroutine mergedj
      integer jjx,jjy,jjz
      forall(jjx=1:nx,jjy=1:ny,jjz=1:nz)
         Dj(jjx,jjy,jjz)=sum(D_j(jjx,jjy,jjz,:))
      end forall
    end subroutine mergedj
  end subroutine proj_sysmat_i2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine sumcs
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    integer::px=0
    select case(isub)
    case(1)
       px=0
    case(2)
       px=1
    case(3)
       px=1
    case(4)
       px=0
    end select
    i=omp_get_thread_num()+1
    call echo_si("SUMCs",i)
    pa=>pas(i)
    pa%ix1=1; pa%iy1=1
    do dy2=oy+0,my-1,2
       iy2=1+dy2
       pa%nsy=my-dy2
       pa%dy=iy2
       do dx2=px+dy2,mx-1,2
          ix2=1+dx2
          sm=>sysmat(ix2,iy2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=ix2
          pa%no=sm%noy
          pa%nsx=mx-dx2
          call cpptra(pa%n,sm%q,i)
          call setma_dr(i)
       end do
    end do
  end subroutine sumcs
  !
  subroutine sumcrs
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMCrs",i)
    pa=>pas(i)
    pa%ix1=1; pa%iy1=1
    do dy2=uy+1,my-1,2
       iy2=1+dy2
       pa%dy=iy2
       pa%nsy=my-dy2
       do dx2=ox+0,dy2-1,2
          ix2=1+dx2
          sm=>sysmat(iy2,ix2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=ix2
          pa%no=sm%noy
          pa%nsx=mx-dx2
          call cpptrb(pa%n,sm%q,i)
          call setma_d(i)
       end do
    end do
  end subroutine sumcrs
    !
  subroutine sumxs
    use omp_lib
    use cuty
    integer i,j
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    integer::px=0
    select case(isub)
    case(1)
       px=0
    case(2)
       px=1
    case(3)
       px=1
    case(4)
       px=0
    end select
    i=omp_get_thread_num()+1
    call echo_si("SUMxs",i)
    pa=>pas(i)
    pa%ix1=mx; pa%iy1=1
    do dy2=oy+0,my-1,2
       iy2=1+dy2
       pa%dy=iy2
       pa%nsy=my-dy2
       if(dy2/=0) then
          j=dy2+px
       else
          j=1+ux
       end if
       do dx2=j,mx-1,2
          ix2=1+dx2
          sm=>sysmat(ix2,iy2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=mx-dx2+1
          pa%no=sm%noy
          pa%nsx=-mx+dx2
          call cpptrax(pa%n,sm%q,i)
          call setma_dr(i)
       end do
    end do
  end subroutine sumxs
    !
  subroutine sumxrs
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMxrs",i)
    pa=>pas(i)
    pa%ix1=mx; pa%iy1=1
    do dy2=oy+2,my-1,2
       iy2=1+dy2
       pa%dy=iy2
       pa%nsy=my-dy2
       do dx2=ux+1,dy2-1,2
          ix2=1+dx2
          sm=>sysmat(iy2,ix2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=mx-dx2+1
          pa%no=sm%noy
          pa%nsx=-mx+dx2
          call cpptrbx(pa%n,sm%q,i)
          call setma_d(i)
       end do
    end do
  end subroutine sumxrs
  !
  subroutine sumys
    use omp_lib
    use cuty
    integer::px=0
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    select case(isub)
    case(1)
       px=0
    case(2)
       px=1
    case(3)
       px=1
    case(4)
       px=0
    end select
    i=omp_get_thread_num()+1
    call echo_si("SUMys",i)
    pa=>pas(i)
    pa%ix1=1; pa%iy1=my
    do dy2=uy+1,my-1,2
       iy2=1+dy2
       pa%dy=my-dy2+1
       pa%nsy=-my+dy2
       do dx2=px+dy2,mx-1,2
          ix2=1+dx2
          sm=>sysmat(ix2,iy2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=ix2
          pa%no=sm%noy
          pa%nsx=mx-dx2
          call cpptray(pa%n,sm%q,i)
          call setma_dr(i)
       end do
    end do
  end subroutine sumys
    !
  subroutine sumyrs
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1
    call echo_si("SUMyrs",i)
    pa=>pas(i)
    pa%ix1=1; pa%iy1=my
    do dy2=uy+1,my-1,2
       iy2=1+dy2
       pa%dy=my-dy2+1
       pa%nsy=-my+dy2
       do dx2=ox+0,dy2-1,2
          ix2=1+dx2
          sm=>sysmat(iy2,ix2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=ix2
          pa%no=sm%noy
          pa%nsx=mx-dx2
          call cpptrby(pa%n,sm%q,i)
          call setma_d(i)
       end do
    end do
  end subroutine sumyrs
    !
  subroutine sumxys
    use omp_lib
    use cuty
    integer i
    integer::px=0
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    select case(isub)
    case(1)
       px=0
    case(2)
       px=1
    case(3)
       px=1
    case(4)
       px=0
    end select
    i=omp_get_thread_num()+1
    call echo_si("SUMxys",i)
    pa=>pas(i)
    pa%ix1=mx; pa%iy1=my
    do dy2=uy+1,my-1,2
       iy2=1+dy2
       pa%dy=my-dy2+1
       pa%nsy=-my+dy2
       do dx2=px+dy2,mx-1,2
          ix2=1+dx2
          sm=>sysmat(ix2,iy2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=mx-dx2+1
          pa%no=sm%noy
          pa%nsx=-mx+dx2
          call cpptraxy(pa%n,sm%q,i)
          call setma_dr(i)
       end do
    end do
  end subroutine sumxys
  !
  subroutine sumxyrs
    use omp_lib
    use cuty
    integer i
    integer ix2,iy2,dx2,dy2
    type(sysmatel),pointer::sm
    type(proj_arg),pointer::pa
    i=omp_get_thread_num()+1      
    call echo_si("SUMxyrs",i)
    pa=>pas(i)
    pa%ix1=mx; pa%iy1=my
    do dy2=oy+2,my-1,2
       iy2=1+dy2
       pa%dy=my-dy2+1
       pa%nsy=-my+dy2
       do dx2=ux+1,dy2-1,2
          ix2=1+dx2
          sm=>sysmat(iy2,ix2)
          pa%n=sm%n
          if(pa%n==0) cycle
          pa%dx=mx-dx2+1
          pa%no=sm%noy
          pa%nsx=-mx+dx2
          call cpptrbxy(pa%n,sm%q,i)
          call setma_d(i)
       end do
    end do
  end subroutine sumxyrs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine init_subset(is)
    integer,intent(in)::is
    isub=is
    if(is==0) return
    ox=xoff(is)
    oy=yoff(is)
    ux=abs(ox-1)
    uy=abs(oy-1)
  end subroutine init_subset
  !
end module projectors
